/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 Mario Orsi
   This file is part of Brahms.
   Brahms is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
   Brahms is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
   You should have received a copy of the GNU General Public License along with Brahms.  If not, see <http://www.gnu.org/licenses/>. */

/*******************************************************************************************************************
 * startUp.c -- system initialization, either by generation 'from scratch' or by loading from previous simulation. *
 *******************************************************************************************************************/

#include "dataStructs.h"

extern VecR *siteCoords, *siteForce, *siteOrientVec, *siteTorque, *siteVelocity, *siteAngMom;
extern RMat *siteRotMat;
extern VecR *atomCoordinates; 
extern Prop boxVolume, xyArea;
extern Site *site;  /* array of Site, represents the state of the system */
extern Atom *atom;  /* array of Atom, represents solute atoms */
extern real *mass, *inertia, *rLipid, timeNow, soluteMass, totalMass;
const int adjustRegion, resetTime;
extern VecR region, solInertia;
extern const real deltaT, extTemperature;
extern RMat boxMatrix;  
extern int nLipids, nTypes, stepCount, loadStructure, loadVelocities, nSymmRigBodies, nPointMasses, nAtoms;
extern const int randSeed, centerInputStruct, nSites, nWaters, nDOPCsDSPCs, nDOPEs, reCenterBilayer, nSolutes;
extern const VecI cells;
extern int *cellList, *interactionType, *nebrTab, *nWatsPerWppSlab;
extern const int diffusion, latDiff, nebrTabMax, rdfWat, maxEdgeCells, applyBarostat;
extern TBuf *tBuf;
extern real *rrDiffuseAv, *ddDiffuseAv, *latSqrDiffAv;
extern const int nValDiffuse, nBuffDiffuse, stepDiffuse, stepWatRdf, limitWatRdf;
extern LatDiffBuff *latDiffBuff;
extern const int nValLatDiff, nBuffLatDiff, stepLatDiff, limitLatDiffAv, zConstraint;
extern const int doCheckpoint, newRun, limitDiffuseAv, writeAreaVol, flexBox;
extern real *histWatRdf, *totalHistEpp, *histPress_xx, *histPress_yy, *histPress_zz;
extern real *histWpp, *soluteHistEpp, *watHistEpp, *headgroupHistEpp, *glycerolHistEpp, *esterHistEpp;
extern int sizeHistWatRdf, countWatRdf;
extern real *soluteHistEdp, *watHistEdp, *cholHistEdp, *phosHistEdp, *glyHistEdp, *estHistEdp, *totalHistEdp;
extern real *CH2_3HistEdp, *CHCH2_2HistEdp, *CH3CH2_2HistEdp;
extern real *totalHistLpp, *histLppTau1, *histLppTau2;
extern const int edp, sizeHistEdp, limitEdp, stepEdp, lpp, sizeHistLpp, limitLpp, stepLpp;
extern const int wpp, sizeHistWpp, limitWpp, stepWpp, epp, sizeHistEpp, limitEpp, stepEpp; 
extern VecR *currWatCoords, *prevWatCoords, *refWatCoords;
extern VecRTwo *xyLipidCom, *xyLipidComTrue;
extern VecRTwo xySoluteCom, xySoluteComTrue;
extern FILE *fPtrZCF, *fPtrAV;
 
/* external function declarations: */
void AccumProps();
void BuildNeighborList();
void BuildQuatFromRotMat( Quat*, RMat* );
void CheckStructureLoaded( int, int );
void CheckStructTypeConsistency( int, int, int, int, int ); 
void ComputeLipLatCom( VecRTwo* ); 
void ComputeTotalMass();
void GetCheckpoint();
void HydrateLipidBilayer( int, int, VecR, real, int );
void InitAngVelsSRB( real, FILE* );
void InitCoords( VecR, FILE* );
void InitDiffusion();
void InitSymmRigBodOrientations( FILE* );
void InitLatDiff(); 
void InitLipidBilayer( int, int, int, VecR, real* );
void InitRand( int );
void InitVels( real, FILE* );
void LoadSoluteParameters( int, int );
void LoadStructureFile( int, int );
void LoadVelFile( int );
void LoadTopology( FILE* );
void PrintBrahmsVersion();
void PrintHeader( FILE* ); 
void PrintPdb( FILE* );
void PrintZConstraintHeader();
void PrintMatrix( RMat );
void PutCheckpoint();
void PutInitialConfiguration();
void RotMatToQuat( Quat*, RMat* );
void SetupCheckpointFiles( FILE* ); 

// global variables:
int DOFs; // total number of degress of freedom of the system's particles - computed in ComputeDOFs();
int countEdp = 0;
int countLpp = 0;
int countWpp = 0;
int countEpp = 0;
real meanEdpSlabVol = 0.;
real meanLppSlabVol = 0.;
real meanEppSlabVol = 0.;
real meanEppSlabArea = 0.;
int halfSizeHistEdp;
int halfSizeHistLpp;
int halfSizeHistWpp;
int halfSizeHistEpp;
real edpDeltaZ; // actual slab thickness for EDP computation
real lppDeltaZ; // actual slab thickness for LPP computation
real wppDeltaZ; // actual slab thickness for WPP computation
real eppDeltaZ; // actual slab thickness for EPP computation
real meanEdpDeltaZ = 0.;
real meanLppDeltaZ = 0.;
real meanWppDeltaZ = 0.;
real meanEppDeltaZ = 0.;
real meanEdpBoxHeight = 0.;
real meanLppBoxHeight = 0.;
real meanWppBoxHeight = 0.;
real meanEppBoxHeight = 0.;
real zConstraintDist; // z-distance at which the constrained solute is located; by convention, z=0 is at the bilayer center
real totMassExceptConstrSolute = 0.;
const VecR regionAdjusted;  // input parameter, used to modify the size of the simulation region */ 

static void AdjustRegion() // Adjusts dimensions of simulation region. 
{                          // The current "region" variable is replaced by "regionAdjusted" (an input parameter from file).
  VSCopy( region, Angstrom_IN_BRAHMS_UNITS, regionAdjusted );
  boxMatrix.u[0]=region.x; boxMatrix.u[4]=region.y; boxMatrix.u[8]=region.z;
  boxMatrix.u[1]=boxMatrix.u[2]=boxMatrix.u[3]=boxMatrix.u[5]=boxMatrix.u[6]=boxMatrix.u[7]=0.;
  boxVolume.val = MatDet( boxMatrix.u );
  xyArea.val = region.x * region.y;
  printf( "Adjusted region:\t   ( %6.3f x %6.3f x %6.3f ) nm\n", region.x * LENGTH_nm, region.y * LENGTH_nm, region.z * LENGTH_nm ); 
}

static void AllocArrays( FILE *fPtr) // Carries out all dynamic array allocations [Rapaport, p.30] 
{ 
  int nb;

  fprintf(fPtr, "AllocArrays... "); fflush(fPtr);

  /* macro AllocMem defined in "defs.h" */
  AllocMem( siteCoords, nSites, VecR );
  AllocMem( siteOrientVec, nSites, VecR );
  AllocMem( siteForce, nSites, VecR );
  AllocMem( siteTorque, nSites, VecR );
  AllocMem( siteVelocity, nSites, VecR );
  AllocMem( siteAngMom, nSites, VecR );
  AllocMem( siteRotMat, nSites, RMat );
  AllocMem( atomCoordinates, nSites, VecR );
  AllocMem( mass, nTypes, real );
  AllocMem( inertia, nTypes, real );
  AllocMem( site, nSites, Site );
  AllocMem( atom, nAtoms, Atom );
  AllocMem( interactionType, nebrTabMax, int );
  AllocMem( nebrTab, 2 * nebrTabMax, int );
  
  if ( nLipids && latDiff ) { // if lipid lateral diffusion calculation is required
    AllocMem( xyLipidCom,     nLipids, VecRTwo ); 
    AllocMem( xyLipidComTrue, nLipids, VecRTwo ); 
  }
  if ( nLipids ) { // allocate arrays to detect transmembrane water permeation - only makes sense for bilayer simulations
    AllocMem( currWatCoords, nWaters, VecR ); 
    AllocMem( prevWatCoords, nWaters, VecR ); 
    AllocMem( refWatCoords, nWaters, VecR );  
  }

  if ( applyBarostat ) AllocMem( cellList, Cube( maxEdgeCells ) + nSites, int );
  else AllocMem( cellList, VProd( cells ) + nSites, int );

  if ( diffusion ) {
    AllocMem( rrDiffuseAv, nValDiffuse, real );
    AllocMem( ddDiffuseAv, nValDiffuse, real );
    AllocMem( tBuf, nBuffDiffuse, TBuf );
    for ( nb = 0; nb < nBuffDiffuse; nb++ ) {
      AllocMem( tBuf[ nb ].orgR, nSites, VecR );
      AllocMem( tBuf[ nb ].rTrue, nSites, VecR );
      AllocMem( tBuf[ nb ].rrDiffuse, nValDiffuse, real );
      AllocMem( tBuf[ nb ].orgD, nSites, VecR );
      AllocMem( tBuf[ nb ].ddDiffuse, nValDiffuse, real );
    }
  }

  if ( nLipids && latDiff ) {
    AllocMem( latSqrDiffAv, nValLatDiff, real );
    AllocMem( latDiffBuff, nBuffLatDiff, LatDiffBuff );
    for ( nb = 0; nb < nBuffLatDiff; nb++ ) {
      AllocMem( latDiffBuff[ nb ].orgTwo, nLipids, VecRTwo );
      AllocMem( latDiffBuff[ nb ].latDispTrueTwo, nLipids, VecRTwo );
      AllocMem( latDiffBuff[ nb ].latSqrDisp, nValLatDiff, real );
    }
  }

  if ( rdfWat ) {
    AllocMem( histWatRdf, sizeHistWatRdf, real );
  }

  if ( edp ) {
    AllocMem( soluteHistEdp, sizeHistEdp, real );
    AllocMem( watHistEdp,  sizeHistEdp, real );
    AllocMem( cholHistEdp, sizeHistEdp, real );
    AllocMem( phosHistEdp, sizeHistEdp, real );
    AllocMem( glyHistEdp,  sizeHistEdp, real );
    AllocMem( estHistEdp,  sizeHistEdp, real ); 
    AllocMem( CH2_3HistEdp, sizeHistEdp, real );
    AllocMem( CHCH2_2HistEdp, sizeHistEdp, real );
    AllocMem( CH3CH2_2HistEdp, sizeHistEdp, real );
    AllocMem( totalHistEdp,sizeHistEdp, real );
    halfSizeHistEdp = sizeHistEdp / 2; // setting variable THIS SHOULD REALLY BE SOMEWHERE ELSE
  }

  if ( lpp ) {
    AllocMem( totalHistLpp, sizeHistLpp, real );
    AllocMem( histLppTau1, sizeHistLpp, real );
    AllocMem( histLppTau2, sizeHistLpp, real );
    AllocMem( histPress_xx, sizeHistLpp, real );
    AllocMem( histPress_yy, sizeHistLpp, real );
    AllocMem( histPress_zz, sizeHistLpp, real );
    halfSizeHistLpp = sizeHistLpp / 2; // setting variable   THIS SHOULD REALLY BE SOMEWHERE ELSE
  }

  if ( wpp ) {
    AllocMem( histWpp,sizeHistWpp, real );
    AllocMem( nWatsPerWppSlab, sizeHistWpp, int );
    halfSizeHistWpp = sizeHistWpp / 2; // setting variable   THIS SHOULD REALLY BE SOMEWHERE ELSE
  }

  if ( epp ) {
    AllocMem( soluteHistEpp,    sizeHistEpp, real );
    AllocMem( watHistEpp,       sizeHistEpp, real );
    AllocMem( headgroupHistEpp, sizeHistEpp, real );
    AllocMem( glycerolHistEpp,  sizeHistEpp, real );
    AllocMem( esterHistEpp,     sizeHistEpp, real );
    AllocMem( totalHistEpp,     sizeHistEpp, real );
    halfSizeHistEpp = sizeHistEpp / 2; // setting variable   THIS SHOULD REALLY BE SOMEWHERE ELSE
  }
  fprintf(fPtr, "done\n"); fflush(fPtr);
}

static void InitForcesTorques() /* Initialize forces and torques */
{
  int n;
  DO_SITE {
    VZero( siteForce[ n ] ); 
    VZero( siteTorque[ n ] ); 
  }
}

static void SetSiteTypes( FILE *fp ) 
{
  int n, nChol = 0, nAmine = 0, nPhos = 0, nEst = 0, nGly = 0, nTail = 0, nWat = 0, nWats, totNSites;

  fprintf(fp, "SetSiteTypes... ");fflush(fp);

  DO_SITE {
    switch ( site[ n ].lipidUnit ) {
    case WATER:     
      site[ n ].type = WATER_TYPE;       
      site[ n ].mechType = SRB; 
      nWat++;  
      break;
    case CHOLINE_OR_AMINE: // this unit can represent both CHOLINE (in DOPC lipids) and AMINE (in DOPE lipids)  
      site[ n ].mechType = PM;   
      if ( site[ n ].lipidSpecies == DOPC ){ 
	site[ n ].type = CHOLINE_TYPE;      
	nChol++; 
      } else if ( site[ n ].lipidSpecies == DOPE ){ 
	site[ n ].type = AMINE_TYPE;      
	nAmine++; 
      }
      break;
    case PHOSPHATE: 
      site[ n ].type = PHOSPHATE_TYPE;      
      site[ n ].mechType = PM;   
      nPhos++; 
      break;
    case GLYCEROL:  
      site[ n ].type = GLYCEROL_TYPE;  
      site[ n ].mechType = SRB;   
      nGly++;  
      break; 
    case ESTER_A:   
    case ESTER_B:   
      site[ n ].type = ESTER_TYPE;   
      site[ n ].mechType = SRB;  
      nEst++;  
      break; 
    case TAIL_A1:  
    case TAIL_A2:  
    case TAIL_A3:  
    case TAIL_A4:  
    case TAIL_A5:  
    case TAIL_B1:  
    case TAIL_B2:  
    case TAIL_B3:  
    case TAIL_B4:  
    case TAIL_B5:  
      site[ n ].type = TAIL_TYPE; 
      site[ n ].mechType = PM; 
      nTail++; 
      break; 
    default: printf( "*** ERROR in SetSiteTypes(), startUp.c: default case.\n" ); 
    }
  }
    
  // DO_SITE printf("## site[ %4d ].lipidUnit = %2d \t site[ %4d ].mechType = %2d\n", n, site[ n ].lipidUnit, n, site[ n ].mechType );
  
  printf( "Structure built: %d waters, %d cholines, %d amines, %d phosphates, %d glycerols, %d esters and %d hydrocarbons\n", 
	  nWat, nChol, nAmine, nPhos, nGly, nEst, nTail );
  totNSites = nWat + nChol + nAmine + nPhos + nGly + nEst + nTail;
  nPointMasses = nChol + nAmine + nPhos + nTail; 
  nSymmRigBodies = nEst + nGly + nWat; 
  nWats = nWat; 
  CheckStructTypeConsistency( totNSites, nSites, nPointMasses, nSymmRigBodies, nWaters );
  fprintf(fp, "done\n");fflush(fp);
}

static void SetPotentialTypes( FILE *fPtr )
{
  int n, ssdCount = 0;
  
  fprintf( fPtr, "SetPotentialTypes... " ); fflush(fPtr);

  DO_SITE {
    site[ n ].mechType = SRB; 
    ssdCount++;
  }

  printf( "There are %d water molecules.\n", ssdCount );
  if ( ssdCount == nSites ) 
    printf( "Simulating 'pure' water system (comprising water sites only)" ); 
  printf("\n");

  fprintf( fPtr, "done\n" ); fflush(fPtr);
}

/************************************************************************************************************
 * CountLipidsPerLayer -- counts the number of lipids in each of the two ("lower" and "upper") monolayers   *
 * making up the whole bilayer. By convention, a lipid molecule belongs to the monolayer where its glycerol *
 * unit is. This function assumes that the origin of the coordinates' system coincides with the centre of   *
 * the simulation region.                                                                                   * 
 ************************************************************************************************************/
static void CountLipidsPerLayer()
{
  int lipidsLowerLayer = 0, lipidsUpperLayer = 0, n;

  DO_SITE {
    if ( GLYCEROL == site[ n ].lipidUnit ) { 
      if ( siteCoords[ n ].z < 0 ) lipidsLowerLayer++;
      if ( siteCoords[ n ].z > 0 ) lipidsUpperLayer++;
    }
  }
  printf("Lipid ditribution: %d molecules in the 'lower' monolayer + %d molecules in the 'upper' monolayer.\n", lipidsLowerLayer, lipidsUpperLayer);
  if ( lipidsLowerLayer != lipidsUpperLayer ) printf("Uneven distribution of lipids in the 'upper/lower' monolayers.\n");
}

/************************************************************************************************************
 * CountWatersPerSection() -- counts the number of water molecules in each of the two halves ("lower" and   *
 * "upper") of the simulation region. This function assumes that the origin of the coordinates' system      *
 * coincides with the centre of the simulation region.                                                      * 
 ************************************************************************************************************/
static void CountWatersPerSection()
{
  int watersLowerSection = 0, watersUpperSection = 0, n;

  DO_SITE {
    if ( WATER_TYPE == site[ n ].lipidUnit ) { 
      if ( siteCoords[ n ].z < 0 ) watersLowerSection++;
      if ( siteCoords[ n ].z > 0 ) watersUpperSection++;
    }
  }
  printf("Water distribution: %3d molecules in the 'lower' region + %3d molecules in the 'upper' region.\n", watersLowerSection, watersUpperSection);
  if ( reCenterBilayer ) printf("\n*** RECENTERING BILAYER ***\n");
}

static void CenterStructCoords()
{
  int n;
  VecR halfRegion;
  
  VSCopy( halfRegion, 0.5, region );
  DO_SITE VVSub( siteCoords[ n ], halfRegion );
}

/************************************************************************************************************
 * ComputeDOFs - computes the number of degrees of freedom in the system                                    *
 ************************************************************************************************************/
static void ComputeDOFs( FILE *fPtr )
{
  int n;

  fprintf( fPtr, "ComputeDOFs... " ); fflush(fPtr);
  
  DOFs = 0; // DOFs is a global variable defined in this module
  DO_SITE {
    DOFs += 3;  // DOFs for a point mass - DOFs as a variable is initialised when defined
    if ( SRB == site[ n ].mechType ) DOFs += 2;  // rotational DOFs of a symmetric rigid body
    else if ( NSRB == site[ n ].mechType ) DOFs += 3;  // rotational DOFs of a non-symmetric rigid body
  }
  DOFs -= 3;   // 3 is the number of additional constraints <- momentum conservation [Rapaport, p.15]
  if ( zConstraint ) DOFs -= 1; // remove 1 DOF as the solute is constrained to a plane z = constant
  fprintf( fPtr, "(Total number of degrees of freedom (DOFs) in the system DOFs = %6d)... ", DOFs );
  fprintf( fPtr, "done\n" ); fflush(fPtr);
}

static void InitWatZArrays() // array employed to detect transmembrane water permeation
{
  int i = 0, n;

  DO_SITE {
    if ( WATER_TYPE == site[ n ].lipidUnit ) { 
      currWatCoords[ i ] = siteCoords[ n ];
      i++;
    }
  }
  DO_WATER refWatCoords[ n ] = currWatCoords[ n ]; 
}

void ZeroEdpHistograms()
{
  int bin;

  for ( bin = 0; bin < sizeHistEdp; bin++ ) { // resetting histograms - to zero
    soluteHistEdp[ bin ] = 0; 
    watHistEdp[ bin ] = 0; 
    cholHistEdp[ bin ] = 0;
    phosHistEdp[ bin ] = 0;
    estHistEdp[ bin ] = 0;
    glyHistEdp[ bin ] = 0;
    CH2_3HistEdp[ bin ] = 0;
    CHCH2_2HistEdp[ bin ] = 0;
    CH3CH2_2HistEdp[ bin ] = 0;
    totalHistEdp[ bin ] = 0;
  }
}   

void ZeroLppHistograms()
{
  int bin;
  for ( bin = 0; bin < sizeHistLpp; bin++ ) { // resetting histograms
    totalHistLpp[ bin ] = histLppTau1[ bin ] = histLppTau2[ bin ] = histPress_xx[ bin ] = histPress_yy[ bin ] = histPress_zz[ bin ] = 0.;
  }   
}

void ZeroWppHistograms()
{
  int bin;
  for ( bin = 0; bin < sizeHistWpp; bin++ ) {
    histWpp[ bin ] = 0.;  // resetting histograms
    nWatsPerWppSlab[ bin ] = 0;
  }   
}

void ZeroEppHistograms()
{
  int bin;
  for ( bin = 0; bin < sizeHistEpp; bin++ ) {
    soluteHistEpp[ bin ]       = 0.;
    watHistEpp[ bin ]          = 0.;
    headgroupHistEpp[ bin ]    = 0.;
    glycerolHistEpp[ bin ]     = 0.;
    esterHistEpp[ bin ]        = 0.;
    totalHistEpp[ bin ]        = 0.;
  }
}

void ZeroWatRdfHistogram() // as in Rapaport p 223
{
  int bin;
  for ( bin = 0; bin < sizeHistWatRdf; bin++ ) histWatRdf[ bin ] = 0;  // resetting histogram
}

static void SetUpDiff()
{
  int stepsPerMeasurement = stepDiffuse * nValDiffuse;
  int stepAvgDiffuse  = stepsPerMeasurement * limitDiffuseAv / nBuffDiffuse;
  real timeWindow = stepsPerMeasurement * deltaT * TIME_ps;
  printf("Diffusion evaluated over %.2f ps measurement times - first average at step %d, then every %d steps.\n", 
	 timeWindow, stepAvgDiffuse + stepsPerMeasurement, stepAvgDiffuse );
  InitDiffusion();
}

static void SetUpLatDiff()
{
  int stepsPerMeasurement = stepLatDiff * nValLatDiff; 
  int stepAvgLatDiff = stepsPerMeasurement * limitLatDiffAv / nBuffLatDiff;
  real timeWindow = stepsPerMeasurement * deltaT;
  printf( "- Lateral diffusion: %.3f-ns windows, first average written to a file \"latDiff-XX.dat\" at %.1f ns, ", timeWindow*TIME_ns, timeWindow*TIME_ns*limitLatDiffAv );
  printf( "then every %.1f ns.\n", stepAvgLatDiff * deltaT * TIME_ns );
  InitLatDiff();
}

static void FindSoluteNAtoms() // open solute.in file and count lines to work out the number of atoms
{
  FILE *fPtr;
  char line[ 81 ];
  int n = 0; // = #lines in the solute.in file
  
  if ( ( fPtr = fopen( "solute.in", "r" ) ) == 0 ) {
    printf( "*** error: Can't open file solute.in!\nExiting...\n" );
    exit( 0 );
  }

  while ( fgets( line, 81, fPtr ) != NULL ) {
    ++n;
  }
  nAtoms = n;
  printf( "The file *solute.in* contains %d lines (atoms).\n", nAtoms );
}

static void InitializeMeasurements( FILE *fPtr )
{
  int anyAnalysis = 0;
  fprintf( fPtr, "InitializeMeasurements... " ); fflush(fPtr);
  AccumProps( 0 );
  
  if ( edp || lpp || wpp || epp || latDiff || diffusion || rdfWat ) anyAnalysis = 1;
  if ( anyAnalysis ) printf("*** The following properties will be measured on-the-fly (as requested via input file \"brahms.an\"): ***\n");
  if ( edp ) { 
    printf( "- Electron densities. Averaged profiles will be written to a file \"edp-XX.dat\" every %.1f ns.\n", stepEdp*limitEdp*deltaT*TIME_ns);
    ZeroEdpHistograms();
    countEdp = 0; /* ??? here or/and somewhere else ??? */
  }
  if ( epp ) {
    printf( "- Electrostatic potentials. Averaged profiles will be written to a file \"epp-XX.dat\" every %.1f ns.\n", stepEpp*limitEpp*deltaT*TIME_ns);
    ZeroEppHistograms();
    countEpp = 0; /* ??? here or/and somewhere else ??? */
  }
  if ( lpp ) { 
    printf( "- Lateral pressure. An averaged profile will be written to a file \"lpp-XX.dat\" every %.1f ns.\n", stepLpp*limitLpp*deltaT*TIME_ns);
    lppDeltaZ = region.z / sizeHistLpp; // computing the binSlab-size i.e. slab thickness - this is used in force.c and lpp.c
    ZeroLppHistograms();
    countLpp = 0; /* ??? here or/and somewhere else ??? */
  }
  if ( wpp ) { 
    printf( "- Water polarization. An averaged profile will be written to a file \"wpp-XX.dat\" every %.1f ns.\n", stepWpp*limitWpp*deltaT*TIME_ns);
    ZeroWppHistograms();
    countWpp = 0; /* ??? here or/and somewhere else ??? */
  }
  if ( latDiff ) SetUpLatDiff(); /* setting up lipid lateral diffusion measurement */
  if ( diffusion ) SetUpDiff(); /* setting up water diffusion measurement */
  if ( rdfWat ) { 
    printf( "- Radial distribution functions. Averaged profiles will be written to a file \"rdf-XX.dat\" every %d steps.\n", stepWatRdf*limitWatRdf );
    ZeroWatRdfHistogram();
    countWatRdf = 0; /* ??? here or/and somewhere else ??? */
    // fPtr = fopen( "rdf.dat", "w" ); fclose( fPtr );
  }
  if ( anyAnalysis ) printf("\n");
  fprintf( fPtr, "done\n" ); fflush(fPtr);
}

void ComputeTotalMass() 
{ 
  int n;
  totalMass = 0.;
  DO_SITE {
    totalMass += mass[ site[ n ].type ]; 
  }
}

void SetupJob( FILE *fp) // Adapted from [Rapaport, pp.20-30]
{
  FILE *fPtr;
  int i;
  real zWaterLayerThickness;
  PrintBrahmsVersion();

  if ( loadStructure ) printf( "Simulation region:  %.2f x %.2f x %.2f \n", region.x, region.y, region.z );

  SetupCheckpointFiles( fp );

  if ( zConstraint || nSolutes ) FindSoluteNAtoms();
  AllocArrays( fp );
  LoadTopology( fp );
  InitializeMeasurements( fp );
  if ( newRun ) {
    if ( loadStructure ) {
      printf("Loading structure from file... \n");
      LoadStructureFile( nSites, nWaters );
      if ( centerInputStruct ) CenterStructCoords();
    } else if ( nLipids || (nWaters<256) ) { // if it's a newRun but there is no structure to load, and if we want a lipid bilayer
      printf("Generating bilayer structure from scratch...\n");
      InitLipidBilayer( nLipids, nDOPCsDSPCs, nDOPEs, region, &zWaterLayerThickness ); // generate bilayer structure from scratch
      HydrateLipidBilayer( nDOPCsDSPCs, nDOPEs, region, zWaterLayerThickness, nWaters ); // hydrate bilayer
      SetSiteTypes( fp );
      InitSymmRigBodOrientations( fp );
    }
  }
  if ( !nLipids ) SetPotentialTypes( fp );
  if ( loadVelocities ) LoadVelFile( nSites );
  ComputeDOFs( fp );
  InitRand( randSeed );
  if ( newRun ) {
    stepCount = 0;
    if ( loadStructure == 0 && stepCount == 0 && nLipids == 0 ) { 
      InitCoords( region, fp ); 
      InitSymmRigBodOrientations( fp );
    }
    if ( loadVelocities == 0 ) {
      InitVels( extTemperature, fp );
      InitAngVelsSRB( extTemperature, fp );
    }
    InitForcesTorques();
    printf (">>>>>>>>>>>>>>>>>>>>> NEW RUN <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n");
    fPtr = fopen( "trajectory.pdb", "w" ); fclose( fPtr );
    if ( doCheckpoint ) { PutCheckpoint(); PutCheckpoint(); }
  } else { // if not new run
    GetCheckpoint();     
    if ( adjustRegion ) AdjustRegion();
    boxVolume.val = MatDet( boxMatrix.u ); // possibly update boxVolume and xyArea variables
    xyArea.val = region.x * region.y;
    ComputeTotalMass();
    if ( resetTime ) { stepCount = 0; timeNow = 0.; }
    CheckStructureLoaded( nSites, nWaters );
    printf( ">>>>>>>>>>>>>>>> CONTINUED RUN <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" );
    printf( "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n" );
  }
  if ( nLipids ) {
    CountLipidsPerLayer();
    CountWatersPerSection();
    InitWatZArrays();
  }
  PrintHeader( stdout );
  //  printf("***sizeHistWpp = %d from SetupJob()\n",sizeHistWpp); fflush(stdout);
  if ( applyBarostat && flexBox && writeAreaVol ) fPtrAV = fopen( "area_volume.dat", "a" ); // opening "global" file, defined on top of this module;

  if ( zConstraint ) {
    for ( i = 0; i < nSites - 1; i++ ) totMassExceptConstrSolute += mass[ site[ i ].type ]; 
    //   printf("totMassExceptConstrSolute = %9.0f amu\n", totMassExceptConstrSolute*MASS_amu); // OK, checked
    fPtrZCF = fopen( "constraintForce.dat", "a" ); // opening "global" file, defined on top of this module;
    PrintZConstraintHeader();
    VSetTwo( xySoluteCom, siteCoords[ nSites-1 ].x, siteCoords[ nSites-1 ].y ); // assumes that the solute is the last element of the "site" array 
    VCopyTwo( xySoluteComTrue, xySoluteCom ); 
  }
  ComputeTotalMass();
  BuildNeighborList();
  PutInitialConfiguration();
}
